!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! the program is coded the process oblique shock wave flow
! properties
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Creator  : B. G.
! Date     : 2015-12-05
! Revised  : 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Note :
!   I. all angle in rad
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine SHOCK(M1_SHOCK,p1_SHOCK,rho1_SHOCK,T1_SHOCK,V1_SHOCK,M1A_SHOCK,eps_SHOCK,&
   &uSBA,vSBA,M2A_SHOCK,thetaS_SHOCK,ps,rhos,Ts,Vs)
implicit none
real*8::M1_SHOCK,p1_SHOCK,rho1_SHOCK,T1_SHOCK,V1_SHOCK,M1A_SHOCK,eps_SHOCK,g
real*8::sinEps,tanEps,M1sinEps2,tanBeta,beta,sinBeta,cosBeta
real*8::uSBA,vSBA,M2A_SHOCK,thetaS_SHOCK,ps,rhos,Ts,Vs

g=1.4
sinEps=sin(eps_SHOCK)                         ! sin(eps)
tanEps=tan(eps_SHOCK)                         ! tan(eps)
M1sinEps2=M1_SHOCK*M1_SHOCK*sinEps*sinEps     ! M1**2*(sin(eps))**2
tanBeta=2*tanEps*(1/M1sinEps2+(g-1)/2)/(g+1)  ! tan(beta)
beta=atan(tanBeta)                            ! beta in rad
sinBeta=sin(beta)
cosBeta=cos(beta)
thetaS_SHOCK=eps_SHOCK-beta
M2A_SHOCK=M1A_SHOCK*sinEps*(2/((g+1)*M1sinEps2)+(g-1)/(g+1))/sinBeta
                                              ! critical Mach number after shock
! non-dimensional velocity  ( like M)
uSBA=M2A_SHOCK*cosBeta      ! critical velocity parallel to shock after shock
vSBA=-M2A_SHOCK*sinBeta     ! critical velocity perpendicular to shock after shock

! flow properties immediately after the shock wave
ps=p1_SHOCK*2*g*(M1sinEps2-(g-1)/(2*g))/(g+1)   ! pressure after shock
rhos=rho1_SHOCK*tanEps/tanBeta                  ! density after shock
Ts=T1_SHOCK*ps*rho1_SHOCK/(p1_SHOCK*rhos)       ! temperature sfter shock
Vs=V1_SHOCK*M2A_SHOCK/M1A_SHOCK                 ! Velocity after shock

end subroutine SHOCK
